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Abstract 

We study both systematic and statistical errors in radiation density ma- 
trix measurements. First we estimate the minimum number of scanning 
phases needed to reduce systematic errors below a fixed threshold. Then, we 
calculate the statistical errors, intrinsic in the procedure that gives the den- 
sity matrix. We present a detailed study of such errors versus the detectors 
quantum efhciency rj and the matrix indexes in the number representation, 
for different radiation states. For unit quantum efhciency, and for both coher- 
ent and squeezed states, the statistical errors of the diagonal matrix elements 
saturate for large n. On the contrary, off-diagonal errors increase with the 
distance from the diagonal. For non unit quantum efhciency the statistical 
errors along the diagonal do not saturate, and increase dramatically versus 
both 1 — 77 and the matrix indexes. 



1 Introduction 

The possibility of investigating quantum radiation states by homodyne de- 
tection techniques recently raised much interest |J . In particular, progress 
has been made on the determination of an exact method to detect the den- 
sity matrix directly from homodyne measurements, in any representation, 
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without resorting to any smoothing or filtering procedure of experimental 
data@, |3], |J. Such method can be basically summarized as follows. By 
means of homodyne detection, the field quadrature = (a^é 1 ^ + ae~ l< ^)/2 is 
measured at any desired phase shift with respect to the local oscillator (a^ 
and a are the creation and annihilation operators of the field mode). Then, 
the density matrix elements are obtained by averaging functions, called "ker- 
nel functions" (or "pattern functions"), on experimental data. We call this 
procedure "homodyning the density matrix" ||, to distinguish it from the 
previously used methods, where the density matrix was reconstructed af- 
ter evaluating the Wigner function as an intermediate step (the celebrated 
"quantum tomography" || |7|, ||). The present method takes into account the 
detectors quantum efficiency, which must be greater than 0.5 for measuring 
the density matrix in the number representation||. 

In this paper, we numerically evaluate the measurement accuracy and 
the statistical errors in homodyning the density matrix, for both unit and 
non unit detectors quantum efficiency rj. In Section || we briefly recali the 
direct method of homodyning the density matrix. Since each matrix element 
is given by an integral over scanning phases, whose number is necessarily 
finite, non negligible systematic errors arise if the number of phases is not 
large enough. Thus, in Section § we numerically estimate a lower value fo for 
the number of phases /, needed for an accurate measurement of a state, and 
we study the dependence of fo on the kind of state. We also investigate the 
convergence of the reconstructed matrix elements towards their respective 
theoretical vàlues as functions of /. In Section [| we introduce the statistical 
errors of the measured matrix elements. We study the errors as functions 
of the matrix indexes and of the quantum efficiency, for both coherent and 
squeezed states. We show that, for r] — 1, the statistical errors of the diagonal 
matrix elements saturate for large n. This result is also analytically obtained 
after introducing an asymptotic approximation for the kernel functions. The 
off-diagonal errors increase with the distance from the diagonal. For 0.5 < 
rj < 1, we show that the statistical errors along the diagonal do not saturate, 
and increase dramatically versus both 1 — rj and the matrix indexes. Due 
to such statistical errors, it is not convenient to use the measured density 
matrix elements to evaluate the expectation vàlues of genèric observables. 
Therefore, in the end of Section [| we consider the possibility of homodyning 
the observable, i.e. measuring directly the expectation value of an observable 
by experiment ally sampling an appropriate kernel function. In particular, we 



2 



consider the number of photons, and we calculate the precision of this kind 
of measurement. In Section |5] we conclude the paper, and in Appendix we 
report some useful calculations in detail. 

2 Homodyning the density matrix 

We briefly recali the method for homodyning the radiation density matrix p. 
Our starting point is the operator identity|| 

= r d0 r+oo ^ Jfe| e _ lkx ^ ? (1) 

JO 7T J-oo 4 

where the field quadrature is defined as x^ = [a^e 1 ^ + ae~ % ^)/2 and is 
the phase with respect to the local oscillator. The trace in Eq.([ï]) can be 
written in terms of quadrature probability distributions p v (x, 0) at phase 0: 
for detectors quantum efficiency rj < 1, such distributions are related to the 
ideal one (r] = 1) by a Gaussian convolution so that in terms of characteristic 
functions one has 

Tr[pe ik ^} = e~ k / dxp v {x, 0) é kx . (2) 

J — oo 

After exchanging integrals over k and x, Eq.([ï|) reads 

P= í— í +00 àxp n {xA)kf{x). (3) 

JO 7T J-oo 

In Equation (|3]), the kernel operator K^\x) is 

Èf{x) = e i ^ a u {ri) {x)e~ i ^ a (4) 

with 

*>&>(*) = f °° dk l 4 e-^ fc2+ ^ : : (5) 

J— oo 4 

(where : : denotes normal ordering). The operator 0^ n \x) can also be written 

as 

0^{x) = d x fi {n \x), (6) 
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with 

fi M {x) = V2 X : : e~ 2 ^ dte* 2 (7) 

Jo 

and x — \fv/ (2?7 — 1). Notice that, equivalently, one has 

^\x) = V2 X e-^e-^e- 2 ^ 2 f^* dte* , (8) 

J o 

where the antidiffusion operator exp(— d 2 /8) is due to normal ordering in 
Eq.©. 

The matrix elements p(n,m) = (n\p\m) are evaluated by averaging the 
kernel functions [i.e. the matrix elements of the kernel operator K^ 1 (x)] 
calculated for random homodyne outcomes. As the experimental data are 
distributed according to the probability p v (x,(p), such average gives a mea- 
surement of the density matrix. In other words, the density matrix is mea- 
sured by experimentally sampling the kernel functions. 

The kernel functions for homodyning the density matrix are reported in 
the following. We have carried out our analysis in the number representat ion, 
for 7] greater than the lower bound 0.5 (it has been shown that rj = 0.5 is a 
universal lower bound for any representation|5|]). 



2.1 Unit quantum efficiency 

For r) — 1, Eq. (|7]) reads: 



£(*) = V2 : e- 3 ^ 9 * : e~ 2 * 2 í 

Jo 



V2x 



dté 2 



(9) 



A simple and fast algorithm can be derived after writing the matrix ele- 
ments (n\fi(x)\m) in factorized form. This technique was first introduced 
by RichterQ for diagonal matrix elements, and was later generalized to off- 
diagonal matrix elements by Leonhardt et aZ.jnj. In Appendix we present a 
simple and alternative derivation which, in our opinion, is useful for further 
developments. The kernel functions are calculated from Eq . (|32D and they 
read 



(m + d\K { f(x)\m) 



Axu m {x)v m+d {x) - 2\fmTïu m+ i{x)v m+d {x) 



2Vm + d+ lu m (x)v m+d+ i\ 



x, 



(10) 
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where the functions Uj(x) and Vj(x) are respectively the normalizable and the 
non normalizable eigenfunctions of the harmònic oscillator (corresponding to 
eigenvalue j). 



2.2 Non unit quantum efficiency 

For i] < 1, no factorization algorithm is known at present. In this case, from 
Eq.(E]) we obtain the following form for the kernel functions[P|: 



(m + d\Kf(x)\m) = eÍ#2 X d+2 y (^)j e " X S (") 



m ( \v 



xJl^iZ-t )i^ + d + l)\x 2u Ke{{-i) d D_ {2v+d+2) {-2i X x)} 



v=0 



where Dj(a) denotes the parabòlic cylinder function. 

3 Systematic errors 

In Equation (|3D the density matrix is given by an integral over the phase <j) 
with respect to the local oscillator. In order to avoid any systematic error, 
one should homodyne the density matrix at perfectly random phases. This is 
the case of the experimental method of Munroe et al.\T3\, where the photon 
number probability distribution is measured by homodyne detection: in such 
measurement no knowledge of the phase is needed, because the diagonal 
kernel functions are independent of <fi. However, for measuring off-diagonal 
matrix elements the knowledge of the random phase is essential, and it is 
difficult to achieve. In such situation, the phase integral is usually performed 
by a phase scanning, as in Ref. f7j. An insufficient number of phases generates 
systematic errors, leading to vàlues for the density matrix elements that are 
far from the true vàlues. Therefore, in the experimental determination of the 
density matrix one has to eliminate these systematic errors as the first step. 

The criterion adopted here to establish the degree of accuracy in a mea- 
surement is based on the absolute deviation of the measured matrix elements 
from the "true" matrix elements. For each p(n,m), obtained from Eq.(|3]), 
we calculate the absolute deviation 

e(n,m) = \p(n,m) - p t (n,m)\ , (12) 
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where pt(n,m) is the true (theoretical) density matrix. For fixed state, the 
set {e(n, m)} (n, m = 0, 1, ...) depends on the number of scanning phases / 
used in the experiment (the number of experimental data per scanning phase 
is kept fixed) . We have an accurate matrix measurement when the maximum 
deviation is reduced below a fixed threshold, for example 

e = maxíeín, m)} < 10~ 4 (n, m — 0, 1...) . (13) 

n,m Lv/J v ' ' \ / 

Now, let us show how the accuracy depends on / for different radiation states. 

The measurement accuracy increases with /. We expect that the more a 
radiation state is either displaced or "asymmetrically" distributed in phase 
space, the higher the number / must be. This is indeed the case. In Fig.[T| we 
show the minimum number of phases fo needed for an accurate measurement 
of coherent and squeezed states. This number increases with both the average 



number of photons (n) and the squeezing parameter r[12 . 

We point out that far off-diagonal kernel functions oscillate very fast as 
functions of 0, thus the larger the matrix dimension, the larger fo. However, 
the main result, i.e. the increase of fo with (n) and r, does not change. 
Indeed, both an increase and a decrease of the matrix dimension merely shift 
the plot in Fig.[I] towards either higher or lower vàlues of fo. In the following 
we set n max = 47. 

A comment about our choice for the accuracy criterion is now in order. 
Our purpose is to show the dependence of f on the average energy and on 
the "asymmetry" in the phase space. This is achieved by calculating the 
absolute deviations {e(n, m)}: indeed, the systematic errors are independent 
of the size of the theoretical matrix element. 

Finally, we briefly examine the dependence on the number of phases / 
for measurements of individual matrix elements. 

We expect that for off-diagonal matrix elements the number of phases 
needed for an accurate measurement is larger than for diagonal ones, due to 
faster oscillations of the integrand in Eq.(§) versus <fi. For coherent states 
this is generally true, as shown for example in Fig.|2], where e(5, 5) < 1CT 4 for 
/ >14, and e(18,5) < 10~ 4 for / >24. For squeezed states the behavior on 
the distance from the diagonal is more complicated. In many cases the same 
result of coherent states is found, see for example Fig.|, where the diagonal 
element p(5, 5) converges faster than p(10,5) and p(15,5) for large enough 
/. However, there are exceptions to this behavior. As an example, in Fig.H 
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we show the asymptotically slower convergence of p(10, 10) with respect to 
p(10,9) and p(10,0). 

4 Statistical errors 

The statistical errors on the measured matrix elements are calculated in 
terms of the errors on real and imaginary parts of the matrix. For a matrix 
element p(n, m) the real part of the statistical variance is defined as 

Re 2 {cx(n, m)} = / — / dxp„(x, 0) Re 2 {(ri|fQ (x)\m)} 

JO 71 J-oo 

- Re 2 {p(n,m)} (14) 

and analogously for the imaginary part. The experimental error of the mea- 
surement is obtained by rescaling the amplitudes \a(n, m) | by a factor 1/yN, 
where N is the total number of experimental data. For simplicity, hereafter 
the quantity a(n,m) will be called statistical error. The statistical errors 
turn out to be independent of / if / > fç,. Thus, we focus attention on the 
general features of the set {cr(n, m)} for different radiation states, at fixed /. 
First we show the results for unit quantum efficiency rj, later we will consider 
the dependence on rj. 

4.1 General features for unit quantum efficiency 

For coherent and squeezed radiation states, the real and imaginary parts of 
the statistical errors exhibit a similar behavior as functions of the matrix 
indexes (with the major exception of the matrix diagonal, where obviously 
Im{cr(n, n)} = 0). Thus, without loss of generality, we can show our results 
in terms of the amplitudes \cr(n, m)\. 

In Fig.|5] we report the matrix of errors \cr(n, m)\ for a coherent state with 
(n) = 4. The contour plot shows that errors increase with the distance d = 
n — m from the diagonal. This is related to the analytical form of the kernel 
operator. In particular, for fixed 0, all the kernel functions are oscillating 
functions of iQ. Moreover, for increasing d the oscillations become faster 
and the oscillation range slowly increases. If a kernel function oscillates fast, 
its statistical average becomes more sensitive to fluctuations of experimental 
data and, therefore, the statistical errors must increase versus d. 
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The contour plot also emphasizes the "saddle region" around the diagonal, 
suggesting that the statistical errors for measured diagonal matrix elements 
saturate to a value independent of n for large enough n. This is shown 
more clearly in Fig|| Such a remarkable feature is general. In fact, it 
is independent on the energy (n) and, more important, it holds for any 
state. Noticeably, the limiting value \a(n, n)\ = \[2 does not depend the 
degree of squeezing. The reason for such saturation is due to the analytic 
form of the diagonal kernel functions. Indeed, the larger n is, the faster the 
kernel functions oscillate vs. x and the errors must increase with n. On the 
other hand, for d = the range of oscillation is fixed between —2 and 2, 
thus the diagonal errors are bounded, and hence they must saturate. These 
considerations are confirmed by considering the explicit form of the statistical 
errors, as given by Eq . (|Ï4D . In particular, from Eq.(^) we can extract the 
relevant contribution for large n upon considering that the kernel functions 
oscillate fast in the region where p(x, 4>) is sizeable. Moreover, p(x, 0) has a 
Gaussian decay, whereas the kernel functions decrease as a power of x. Thus, 

í r d(f) f + °° 1 
\a(n, n)\ ~ < J — J dxp(x, <p) 4 cos 2 (k n x) > . (15) 

For large vàlues of n, k„ — > oo: if p(x, <p) can be considered constant over a 
cycle Ax = n/k n , the integral over x in Eq. ( [15] ) gives just the average of 
cos 2 (k n x), which leads to 

\(r(n,n)\ ~ V2. (16) 

If very squeezed states are considered [i.e. with very sharp p(x, <$)] the errors 
will saturate for larger n. In Fig.^we show \a(n, m) \ for a squeezed state with 
(n) = 4 and r = 1: the plot is quite different from Fig.[|, but the diagonal 
errors still saturate. 

4.2 Dependence on the quantum efficiency 

The influence of the quantum efficiency rj on \a(n,m)\ is very strong. In- 
deed, for non unit quantum efficiency of detectors the behavior of the kernel 
functions (|Ï2|) changes dramatically: for fixed n and m, the oscillation range 
increases very rapidly as r] approaches the lower bound r] = 0.5, and the 
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resulting errors increase rapidly as well. The growth rate is different for dif- 
ferent matrix elements: as an example, in Fig.[|]we show some diagonal errors 
as functions of quantum emciency (for a coherent state). Furthermore, the 
diagonal errors no longer saturate for large vàlues of n. Very similar results 
are found for squeezed states. In particular, the growth rate of diagonal er- 
rors \a(n, n) | vs. 1 —77 is slightly larger than for coherent states. The diagonal 
errors for a squeezed state are shown in Fig|| for different vàlues of 77. 

For fixed 77 < 1, the oscillation range of the kernel functions increases 
with both n and the distance d from the diagonal. Thus, for increasing n 
and d the statistical errors increase. For example, we consider 77 = 0.99: after 
a comparison between Fig.|H] and Fig.[5], one can see that the open contour 



levels for \cr(n, m)\ close, and any errors saturation disappears. Figure |TÜ 



shows that dràstic modifications arise with respect to the ideal case for 77 = 1. 
This means that, already for 77 = 0.99, in order to have the same experimental 
errors on the measurement of the density matrix, the number of data must 
be much larger than in the ideal case. 



4.3 Precision of homodyning observables 

From the measured density matrix, one can evaluate the probability distri- 
butions of operators that are functions of the field operators a and a) . Thus, 
by means of homodyne experimental data it is possible to obtain indirect 
measurements of observables. However, for some observables the propaga- 
tion law of statistical errors leads to additional noise with respect to direct 
detection. In some cases, such indirect detection through the density matrix 
can be overcome by a more convenient procedure, namely homodyning the 
observable. By homodyning the observable we mean the experimental sam- 
pling of an appropriate kernel function, which directly gives the expectation 
value of the desired observable. 

We consider, as an example, the homodyne measurement of the mean 
photon number (n). From Eq. ([]]), (n) is expressed as 

(n) = /"^ [ + °° dxp(x,<P)F(x,4)) (17) 

JO 71 J-00 

where 

F(x, 0) = F(x) = f> £J dk M e-t fc2 +^ (Ç\ =2x 2 -\. (18) 
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In Equation (0), L$ denote zero-order Laguerre polynomials and unit de- 
tectors efficiency has been considered. The statistical fluctuations of the 
measured mean photon number are given by 

4> = r — r°° *) p2 ^ <t>) - • ( i9 ) 

w JO TT J-oo 

and c/ft) is the statistical error for homodyning the mean photon number. 
The precision of this homodyne measurement is defined by the relation 

e\ h) = 4> - (An 2 ) , (20) 

where (An 2 ) is the intrinsic quantum uncertainty 

(An 2 ) = (h 2 ) - (n) 2 = (a f V> + (fi) - (n) 2 . (21) 



The uncertainty (An 2 ) can be expressed in terms of quadrature probability 
distributions: after calculating the kernel function for operator a^ 2 a 2 [PH , 
Eq.(gg) reads 

(An 2 ) = f ^ Í +C °dx p(x, <P) (-a; 4 - 2x 2 \ - (h) 2 . (22) 

JO TT J-oo 13 J 

In conclusion, the precision for homodyning the photon number is 

e {h) = ^=((An 2 ) + (n) 2 + (n) + l) 1/2 . (23) 



5 Conclusions 

We analyzed both systematic and statistical errors for homodyne detection 
of the density matrix of light. Such detection is performed by suitably pro- 
cessing homodyne experimental data. We studied the behavior of systematic 
errors as functions of the number of scanning phases /. We calculated the 
lower bound for /, needed for an accurate matrix measurement. We found 
that this lower bound increases with both the mean photon number and the 
"asymmetry" in phase space of the state. Then we considered the statisti- 
cal errors corresponding to the data average that gives each matrix element. 
Noticeably, for unit quantum efficiency of the detectors the diagonal errors 
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<r(n, n) of the matrix elements p(n, n) "saturate" to the fixed value \[2 for 
large enough n. This feature is independent of the degree of squeezing. The 
off-diagonal errors increase with the distance from the diagonal. If detectors 
quantum emciency is decreased, the errors increase for each matrix element. 
In particular, any saturation effects disappear. Finally, we considered the ho- 
modyne detection of the mean photon number, that is achieved by sampling 
an appropriate kernel function, and we analytically evaluated the precision 
of such measurement. We think that the results presented here are relevant 
from a fundamental point of view and provide the experimentalist with im- 
portant information on the behavior of errors in homodyning the density 
matrix. 



The factorization of the matrix element (m + d\p,(x)\m) is performed in two 
steps. By setting n = m + d we obtain 



Then, the derivatives with respect to x and the summation are evaluated as 
follows. We introduce the "seed functions" 



that generate two sets of functions {uj(x)} and {vj(x)} for j = 0, 1,2..., as 



6 Appendix 




(24) 




(25) 



(26) 




(27) 



(28) 
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By means of the following identity between operators 

d x uo(x) = u (x)(d x - 2x) , 



we obtain 



-y) Uq(x)v (x) 



d\ u (x)vd(x) . 



(29) 



(30) 



As noticed in Ref. pH| , the functions {uj(x)} and {vj(x)} are respectively 
the normalizable and the non normalizable eigenfunctions of the harmònic 
oscillator (corresponding to eigenvalue j). Thus, by using the Standard re- 
cursion relations for the harmònic oscillator eigenfunctions, we can easily 
demonstrate the following identity pi 1: 



1 



d x 



2i> 



~ü\\~2) Uo[ - x > Vd[ - x > 



j=0 \ 



J + d 
3 



\v-3 



u + d 
J+d 



Uj(x)v j+d (x) .(31) 



After substituting (p0| ) and (|3l|) in Eq . (pi|) , we obtain the factorized formula 

(m + d\fi(x)\m) = (m\p,(x)\m + d) = u m (x)v m+d {x) , (32) 

where we use the fact that fi(x) is real selfadjoint. Finally, the kernel func- 
tions are obtained using Eqs. and (|32|). 
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Figure 1: Minimum number of scanning phases /o required by the condition 
e < 10~ 4 vs. the mean number of photons (n) for coherent states (circles), 
squeezed states with r = 0.6 (triangles), and r = 1 (squares). [The matrix 
dimensions are fixed to n max = 47.] 

Figure 2: Absolute deviation e(n,m) vs. / for a coherent state with (n) = 
4: (n, to) = (5,5) (circles), (n,m) = (10,5) (triangles), (n, to) = (18,5) 
(squares). The theoretical matrix elements are p t (5, 5) = 0.15629, p t (10, 5) = 
0.02876, #(18,5) = 0.00017. 

Figure 3: Absolute deviation e(n, to) vs. / for a squeezed state with (n) = 4, 
r = 1: (n, to) = (5,5) (circles), (n, to) = (10,5) (triangles), (n, to) = (15,5) 
(squares). The theoretical matrix elements are p*(5, 5) = 0.04182, pt(10, 5) = 
0.03231, #(15,5) = 0.01852. 

Figure 4: Absolute deviation e(n, to) vs. / for a squeezed state with (n) = 4, 
r = 1: (n, to) = (10,10) (circles), (n, to) = (10,9) (triangles), (n, to) = 
(10,0) (squares). The theoretical matrix elements are p f (10, 10) = 0.02495, 
Pí (10,9) = 0.02418, p t (10,0) = 0.09307. 

Figure 5: Statistical error amplitudes \cr(n, m)\ for a coherent state with 
(n) = 4 (77 = 1). 

Figure 6: \cr(n, n)\ for: coherent state with (n) = 4 (circles), squeezed state 
with (n) = 4, r = 1 (triangles) (rj = 1). 

Figure 7: \a(n,m)\ for a squeezed state with (n) = 4, r = 1 (77 = 1). 

Figure 8: \cr(n, n)\ vs. 1 — i] for n = 0,2,5,15 on a semilogarithmic 
scale (for a coherent state with (n) = 4). The quantum eficiències are 
rj = 1,0.99,0.97,0.95,0.9. 
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Figure 9: \cr(n, n)\ for a squeezed state with (n) — 4, r — 1 for: rj = 1 
(circles), ?7 = 0.99 (triangles), = 0.97 (squares), rj = 0.95 (rhombi), rj = 0.9 
(stars). 



Figure 10: \cr(n, m)\ for a coherent state with (n) = 4 for quantum efficiency 
7] = 0.99. 
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Fig.l G.M. D'Ariano et a 
Systematic and statistical errors.. 




Fig.2 G.M. D'Ariano et al 
Systematic and statistical errors.. 




Fig.3 G.M. D'Ariano et al 
Systematic and statistical errors.. 




Fig.4 G.M. D'Ariano et a 
Systematic and statistical errors.. 




Fig.5 G.M. D'Ariano et a 
Systematic and statistical errors.. 
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Fig.8 G.M. D'Ariano et al 
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Fig.9 G.M. D'Ariano et a 
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